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Abstract 

Crystalline bilayers of charged colloidal suspensions which are confined between two parallel 
plates and sheared via a relative motion of the two plates are studied by extensive Brownian dy- 
namics computer simulations. The charge-stabilized suspension is modeled by a Yukawa pair po- 
tential. The unsheared equilibrium configuration are two crystalline layers with a nested quadratic 
in-plane structure. For increasing shear rates 7, we find the following steady states: first, there is 
a static solid which is elastically sheared until a yield-stress limit is reached. Then there are two 
crystalline layers sliding on top of each other with a registration procedure. Higher shear rates melt 
the crystalline bilayers and even higher shear rates lead to a reentrant solid stratified in the shear 
direction. This qualitative scenario is similar to that found in previous bulk simulations. We have 
then studied the relaxation of the sheared steady state back to equilibrium after an instantaneous 
cessation of shear and found a nonmonotonic behavior of the typical relaxation time as a function 
of the shear rate 7. In particular, application of high shear rates accelerates the relaxation back to 
equilibrium since shear-ordering facilitates the growth of the equilibrium crystal. This mechanism 
can be used to grow defect-free colloidal crystals from strongly sheared suspensions. Our theoretical 
predictions can be verified in real-space experiments of strongly confined charged suspensions. 

PACS numbers: 82.70Dd, 83.10.Mj, 61.20.Ja 
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I. INTRODUCTION 



A fundamental understanding of the different processes governing the relaxation of 
metastable phases back to equilibrium is critical for many basic questions in condensed 
matter physics and material science. Also, relaxational processes are omnipresent in in- 
dustrial applications. Colloidal suspensions represent excellent model systems where such 
questions can be studied directly in real-space as the length scales are conveniently accessed 
experimentally, the (variable) interactions can be described theoretically in a simple way 
and the microscopic processes are rather slow as compared to molecular materials. This 
has extensively exploited in previous studies of interaction dependent equilibrium properties 



and dynamics 



QflflQ. 



One important example for a nonequilibrium steady state is a 
sheared colloidal suspension. It is known that application of shear may destroy the under- 
lying equilibrium crystalline structure of the unsheared suspension and can also lead to 
a reentrance ordering for high shear rates fG] . After cessation of shear the system will relax 
back to equilibrium from the sheared steady state. The microscopic details of this relaxation 
process are far from being resolved. 



If an additional confinement between two parallel plates is considered pj, various ex- 
periments sl, l| 0, 0, 0, 13 reveal a rich and subtle influence of shear on the structure 
close to the wall. Accordingly the relaxation back to equilibrium after cessation of shear is 
a fascinating but complex process which is a competition between wetting effects near the 
walls and bulk relaxation. In experiments on strongly confined suspensions, for instance, 
a complex pathway of the relaxation back to equilibrium was obtained 14j : a bilayer bcc 
crystal was shear-molten to recrystallize as a buckled single layer triangular lattice which 
subsequently underwent a martensitic transition back to the equilibrium phase. 

Most of the theoretical studies on colloidal suspensions have addressed the influence of 
linear shear flow on the bulk structure via non-equilibrium Brownian dynamics (NEBD) com- 
puter simulations 01 where hydrodynamic interaction s are neglected an d involve charged 
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mm 



24. 



colloidal particles modeled by a Yukawa pair interaction 0,0 0,0,0, 

Shear-induced melting of colloidal bulk crystals and subsequent reentrant ordering at higher 
shear rates are confirmed by simulation. However, simulational work including a wall act- 
ing on a sheared suspension is sparse; apart from a NEBD simulation in a channel 
theoretical investigations were for a single colloidal particle only (22], 0|. 
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In the present paper we address the relaxation of shear-induced structures after cessation 
of shear. We use the standard Yukawa model for confined systems and employ NEBD 
simulations. Here we focus on the simple and transparent situation of a colloidal bilayers 
which are confined between two parallel plates and sheared via a relative motion of the two 
plates. The reasons to do so is three-fold: First, the equilibrium phase diagram for confined 
crystalline bilayers interacting via a Yukawa pair potential is known from recent lattice- 
sum techniques at low temperatures j^J. This phase diagram was recently confirmed in 
experiments on charged suspensions strongly confined between two glass plates 3(3] . Second, 
the structure and the defects in a crystalline bilayer are easier to classify than in a multilayer. 
Last not least there are experimental studies for strongly confined situations which are 



not completely understood and are a c hal 
simulation studies of Das and coworkers 



enge for a theoretical treatment Recent 



321 ] have addressed similar questions regarding 



sliding bilayers. The model employed in the studies of Das et al, however, is simpler than 
ours, it does not possess a spatial dimension z perpendicular to the plates and hopping 
processes between the layers ar e ig nored. Furthermore, the relaxation back to equilibrium 
is not investigated in Refs. j^l], Isj ]. 

In order to be specific, we chose the unsheared equilibrium configuration to be two crys- 
talline layers with a nested quadratic in-plane structure. This is the same staring configu- 
ration as used in the experiments |l4j. For increasing shear rates 7, we find the following 
scenario of steady states: first, there is a static solid which is elastically sheared until a yield- 
stress limit is reached. Then there are two crystalline layers sliding on top of each other with 
a lock-in registration procedure similar to that observed in recent experiments by Palberg 
and Biehl p^L [mJ. Higher shear rates melt the crystalline bilayers and even higher shear 
rates lead to a reentrant solid stratified in the shear direction. This qualitative scenario is 
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2l\ . The shear-induced ordering 



similar to that found in previous bulk simulations 
at high shear rates is reminiscent to the transition towards lane formation in oppositely 
driven p artlC ,e S Q. We have then s tn dl ed the Nation of the sheared S tea dy state oacK 
to equilibrium after an instantaneous cessation of shear and found a nonmonotonic behavior 
of the typical relaxation time as a function of the shear rate 7. In particular, application 
of high shear rates accelerates the relaxation back to equilibrium via shear ordering in the 
steady state. This mechanism can be used to grow defect-free colloidal crystals from strongly 
sheared suspensions as was proposed by Clark and coworkers Our theoretical pre- 
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dictions can be verified in experiments of confined charged suspensions 

The paper is organized as follows: In section |Hj we introduce the ground state model 
model for crystalline bilayers. The nonequilibrium Brownian dynamics simulation technique 
is explained in section 11111 Results are presented in section I1VI Finally we conclude in 
section |V| 

II. THE MODEL 

In this part, we define our model. This is basically a generalization towards finite temper- 
ature of the ground state model used in Ref. j^j concerning the equilibrium (i.e., without 
external applied shear flow) phase diagram of crystalline colloidal bilayers interacting via a 
Yukawa potential. In detail, our system consists of two layers containing in total N particles 
in the (x,y) plane. The total area density of the two layers is p = N/A with A denoting 
the layer area in the (x, y) plane. The distance D between the layers in the z direction is 
prescribed by the external potential confining the system. The particles are interacting via 
the Yukawa pair potential 

exp(-M-) 

V yU k{r) = V , (1) 

KT 

where r is the center-center separation. The inverse screening length k which governs the 
range of the interaction is given in terms of the micro-ion concentration by Debye-Hiickel 
screening theory. The energy amplitude Vq = Z 2 exp(2/«i2)«;/e(l + kR) 2 scales with the 
square of the charges Z of the particles of physical hard core radius R reduced by the 
dielectric e permittivity of the solvent (e = 1 for the dusty plasma). Typically, Z is of the 
order of 100 — 100000 elementary charges such that V yu k(r) at typical interparticle distances 
can be much larger than the thermal energy k^T at room temperature justifying formally 
zero-temperature calculations. Then the energy scale is set by Vo alone and phase transitions 
in large bilayer systems are completely determined by two dimensionless parameters, namely 
the reduced layer density: 



V = pD 2 /2 (2) 



and the reduced screening strength: 



A = kD (3) 
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For zero temperature, the stable state is solid but different crystalline structures of the 
bilayers are conceivable. The result for the phase diagram in a (77, A)-map can be found in 
Ref. j^J. Here, we explore the same model for finite temperature by computer simulation. 



III. THE NONEQUILIBRIUM BROWNIAN DYNAMICS COMPUTER SIMULA- 
TION 

A. Simulation method 



Here, we provide a detailed description of our Brownian dynamics method that was used 
to investigate non- equilibrium sheared colloidal bilayers (at finite temperature). A schematic 
setup of the system in the (x, z) plane is depicted in Fig. [T] The integration scheme for our 
model system in the presence of an external steady shear rate 7 reads: 



ii(t + St) = Ti (t) + 7%Fi(t)<Jt + SWi + ^ Zl {t)5te a 
k B T 



(4) 



Thereby r^t) = [xi(t),yi(t),Zi(t)] is the position of the i—th colloidal particle at time t and 
D denotes its free diffusion constant. All the contributions to the equation of motion (J3J) 
are explained below. 

Within a small time interval St, that particle moves under the influence of the sum of 
conservative forces Fj(i) stemming from: (i) The pair interaction V yu \. [see Eq. ([Q] between 
particle i and the neighboring ones, (ii) The repulsive interaction with the soft wall(s) whose 
potential of interaction, V wa u, is modeled as follows 

10 



Vwall{ z ) 



ae LJ 



+ e L j, for 



> 



1/6 



for 



— \z\ 



< 



1/6 



(5) 



where 



a 



10 



3.07002 



1 /6 

[with z min = (|) a L j minimizing V wa u in Eq. ©] so that V wa u(a LJ ) = e L j ■ This 
(truncated and shifted) 10 — 4 Lennard- Jones potential given by Eq. © assumes that we 
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FIG. 1: View in the (x,z) plane of the setup of the colloidal bilayer confined between two soft 
walls. 

have thin soft walls . Note that the use of a 9 — 3 Lennard- Jones potential corresponding to 
semi-infinite walls would not change qualitatively the results. Also the use of charged hard 
walls would not affect our main results. 

Furthermore due to the presence of the solvent, the particles experience (i) a friction 
whose constant is given by /cbT '/ D and (ii) random displacements, 5Wj. Those latter are 
sampled from a Gaussian distribution with zero mean and variance 2Do5t (for each Cartesian 
component). The last term in Eq. (0J) represents the applied shear in the x— direction, and 
imposes an explicit linear flow field. The zero velocity plane of the imposed shear lies at the 
midplane between the plates. 

B. Parameters 

The colloidal particles are confined in a rectangular L x L x D wa u box where periodic 
boundary conditions are applied in the (x, y) directions. The system is made up of N = 800 
particles (i. e., 400 particles per layer). The units are set as follows: ksT = 1/(3 sets the 
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FIG. 2: Intralayer (x,y) pair distribution function g{r = \j x 2 + y 2 ) at equilibrium (7 = 0). The 
inset shows a simulation snapshot where the filled (open) circles represent particles belonging to 
the upper (lower) layer. 

energy scale, the (typical average) interlayer separation D = D wa u — 2au (see also Fig. 

sets the length scale, and r = D 2 /D sets the time scale. For the Yukawa interparticle 
interaction [see Eq. we choose (3Vq = 6000, whereas for the wall-particle interaction [see 
Eq. we choose jSeij = 1 and a^j = O.LD. The time step was set to 5t = 10~ 5 r. The 
reduced colloidal particle density is set to 77 = = 0.24 (so that L = 40.82D) and the 
reduced screening is A = nD = 2.5. Those latter parameters lead to the staggered square 
phase in the ground state (or at very low temperature) as can be seen on the phase diagram 
from Ref. 29|. 

The equilibrium (i.e., 7 = 0) properties of our model system are obtained over a period of 
10 6 BD time steps (i. e., lOr.) The corresponding in-plane (x,y) pair distribution function 
g(r) is shown in Fig. It clearly shows a high degree of ordering as characterized by the 
pronounced peaks and the deep minima. The snapshot also provided in Fig. |2] confirms the 
square lattice structure expected for those parameters. To quantify the layer extension in 
the ^-direction we have also plotted the the particle density n(z) that can be found on Fig. 
01 The mean interlayer separation is then given by 2 j® wal1 ' 2 zn(z)L 2 dz m 0.99D, so that 
(in practice) D corresponds indeed to the interlayer separation. 
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FIG. 3: Laterally averaged inhomogeneous particle density n(z) at equilibrium (7 = 0). 
IV. RESULTS 

A. Effect of shear flow 

Starting from the equilibrium configuration described in the previous section, an external 
shear is applied during a period of 4 x 10 6 BD steps (i.e., 40r). A steady is reached after 
typically lOr, and subsequent measurements are performed over a typical period of 20 r. 

It is instructive to start our study by analyzing the microstructures reported in Fig. H] 
corresponding to different 7. From a structural point of view one can (qualitatively) identify 
three regimes: 

• At sufficiently low shear rates (here 7 = 20/r and 7 = 50/r), it can be seen that the 
crystalline structure (namely square) as well as the degree of ordering are conserved 
compared to the equilibrium situation (i.e., 7 = 0- see Fig. EJ). Consequently we are 
in an elastic regime where the applied shear flow is smaller than the yield stress. 

• For intermediate shear rates (here 7 = 60/r and 7 = 80/r), there is a (relative) strong 
disorder and the structure can therefore be qualified as liquid. In other words we have 
to deal with a shear induced melting. 

• At high shear rates (here 7 = 100/r and 7 = 200/r), the system gets again ordered 
(especially for the highest shear rate 7 = 200/r) but exhibits a different (intralayer) 
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FIG. 4: Simulation snapshots for different values of the shear rate 7 (as indicated) where the filled 
(open) circles represent particles belonging to the upper (lower) layer. 
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FIG. 5: {a) Interlayer (x,y) pair distribution function g(r = \/ x 2 + y 2 ) for small values of 7 
(as indicated in the legend), (b) Intralayer (x,y) pair distribution function g(r = \J x 2 + y 2 ) for 
different values of 7 (as indicated in the legend). The corresponding simulation snapshots are 
displayed in Fig. El 



crystalline symmetry (namely a triangular lattice) than the equilibrium one. Conse- 
quently, we have a reentrant behavior concerning the intralay ^r-ordering upon shear- 
ing. 

In order to obtain a more quantitative description of these 7-dependent structural prop- 
erties, we have also computed the (azimuthally averaged) interlayer- and intralay er-pa.ii- 
distribution-functions g(r = \J x 2 + y 2 ) for different 7. The results are presented in Fig. 
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5. 

The elastic behavior and the yield stress can be best understood by considering the 
interlayer and intralayer g(r). From Fig. we see that at weak shearing (here 7 = 20/ t), the 
intralayer crystalline structure as well as the interlayer-lattice-correlation remain unchanged 
compared to the 7 = case (the latter is not reported in Fig. |5J). At larger shear rate (here 
7 = 50/r) the degree of interlayer-lattice-correlation gets weaker than that of the intralayer 
one. A closer look at Fig. EJa) reveals that, for 7 = 50/r, the first peak is (asymmetrically) 
splitted into two neighboring peaks. This is the signature of a small relative displacement of 
the two square layer-lattices. Upon further increasing the shear rate (now at 7 = 60/r), the 
bilayer becomes a liquid, demonstrating that there is a yield stress 70 whose value is such 
that 50/r < 70 < 60/r. 

Above the yield stress, the intralayer g(r) [see Fig. GJb)] exhibits a non-trivial behavior 
with respect to the applied shear flow, in agreement with our previous discussion on the mi- 
crostructures depicted in Fig. HJ More precisely, at intermediate values of 7 (here 60/r and 
80/r), the intralayer layer structure corresponds to a liquid one. Nonetheless and interest- 
ingly, at first neighbor separations, the square structure locally persists, but in coexistence 
with a triangular structure, as indicated by the broadened (splitted) first peak. This feature 
can also be nicely visualized on the snapshots from Fig. HJ At high shear rates (here 100/r 
and 200/r), there is a strong (re)ordering into a triangular lattice as indicated by the shifted 
first pronounced peak (especially for 7 = 200/r). However the degree of ordering reported 
for those highly sheared structures is not as high as that observed below yield stress. 

To further quantify the behavior of highly sheared colloidal bilayers and also to provide 
a dynamical information, we are going to examine the (dimensionless) modified Lindemann 
parameter, Ti,(t), that is defined as follows 

rut) - <if>, (6) 

where (u 2 (t)) corresponds to the difference in the mean square displacement of neighboring 
particles from their initial sites r = r(t = t ). More explicitly, (u 2 (t)) can be written as 



(u 2 (t)) 



> N N b \ 

n E w E [(*(*) - r ^°)) - ^ - r ^°))] 2 ) > ( ? ) 



where Ti(t) = [xi{t),yi(t)], (. . . ) denotes an averaging over BD steps, the index j stands for 
the Nt, nearest neighbors of particle i lying in the same upper or lower layers. Typically, for 
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FIG. 6: Modified Lindemann parameter Ti(t) for different values of 7 (as reported in the legend) 

a (local) triangular lattice environment Nf, = 6 while for a rectangular one iV& = 4. Besides, 
we also average over several reference times to to improve the statistics. Due to the finite 
size of the simulation box, one is typically limited to observation times At b s of the order of 
At obs « Lj^ max D « 0.2r (by taking here ^ max = 200/r). 

Our results are presented in Fig. El In the elastic regime (small 7), the Lindemann param- 
eter T(t) exhibits a plateau at "long" times confirming the crystalline intralayer structure. 
At higher 7 (i.e., 7 > 60) the situation gets more complicated. For 60/r < 7 < 100/r, T L (t) 
diverges proving a liquid behavior. While this feature was clearly expected for 7 = 60/r and 
80/r from our static analysis of g(r) [see Fig. Efb)], that was not obvious for 7 = 100/r. 
It is therefore only at very high shear rate (i.e., 7 > 200) that true intralayer crystalline- 
reordering is recovered, as indicated by the existence of the plateau in Ti(t) whose value is 
comparable to that obtained in the elastic regime. 

B. Relaxation after cessation of shear 

We now investigate how the system gets back to equilibrium after cessation of shear. A 
suitable and simple way to study a relaxation process is to monitor the evolution in time 
of the total potential energy of interaction E(t) = V yu k + V wa u. In our simulations, the 
cessation of shear occurs at t = 40r. Profiles of E(t) for different shear rates 7 applied prior 
relaxation are plotted in Fig. The corresponding microstructures at long time t = 80r 
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FIG. 7: Time evolution of the total potential energy of interaction E(t): before - during - and after 
shear. The values of 7, considered during the shear process, are reported in the legend. 

for 6O7/T < 7 < 2007/r are sketched in Fig. |H| For low 7 (here 7 < 50j/t), the relaxation 
process is very fast as it should be. Note that the equilibrium energy value is not exactly 
recovered because of the existence of some long-living defects. 

The relaxation process gets qualitatively different for highly sheared systems (here 7 > 
607/r). For the samples that have undergone a shear- induced melting as deduced from our 
criterion based on Ti{t) [see Fig. El with 7T = 60, 80, 100], we remark that they all exhibit a 
similar relaxation behavior [see Fig. 0with 7T = 60, 80, 100]. In particular the relaxation is 
thereby much slower, partly due to the existence of many long living defects. Those latter 
also explain the high energy reported in the long time scale. There are several defects such 
as dislocations, (low angle) grain boundaries (especially for 7T = 60, 80) and vacancies that 
are easily identifiable in the snapshots of Fig. |H1 

On the other hand at large enough 7 (here 7 = 200/r) the relaxation is faster as indicated 
by the faster earlier occurrence of a .£?(£) -plateau (which is also deeper). Nonetheless, the 
energy of this (nearly) relaxed system remains higher than those that were weakly sheared 
(7 < 7o)- Again the existence of some vacancies (see Fig. |H]with 7 = 200/r) increases the 
energy system as well as the time of full relaxation. 
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FIG. 8: Simulation snapshots of relaxed systems taken at t = 80r for different values of the prior 
applied shear rates 7 (as indicated). 

V. CONCLUSIONS 

To conclude we perform Brownian dynamics computer simulations to study crystalline 
bilayers of charged colloidal suspensions which are confined between two parallel plates and 
sheared via a relative motion of the two plates. For the parameters under consideration, 
the unsheared equilibrium configuration are two crystalline layers with a nested quadratic 
in-plane structure. For increasing shear rates 7, we find the following steady states: first, 
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there is a static solid which is elastically sheared until a yield-stress limit is reached. Higher 
shear rates melt the crystalline bilayers and even higher shear rates lead to a reentrant 
solid stratified in the shear direction. We have then studied the relaxation of the sheared 
steady state back to equilibrium after an instantaneous cessation of shear and found a 
nonmonotonic behavior of the typical relaxation time as a function of the shear rate 7. In 
particular, application of (very) high shear rates accelerates the (post-)relaxation back to 
equilibrium since shear-ordering facilitates the growth of the equilibrium crystal. 
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